disp('hqlyzdl');

a = [1,2];
b = [2,1];
psi = fn_lie_oneterm(a,b);
disp(psi.amp);
disp(psi.base);

Num = 80;
Numb = 40;
bq = cell(1, Numb);
bdq = cell(1, Numb);

for i = 1:Numb
    bq{i} = fn_bq(i,Num);
    bdq{i} = fn_bdq(i,Num);
end


psi = fn_lie(bq{3},bdq{3});
fn_print(psi);

fn_plot(psi,1);

for k = 1:5
    psi = fn_lie(psi,bdq{3});
    fprintf('k = %d \n', k);
    fn_print(psi);
end

fn_print(psi);

fn_plot(psi, 4);

tpsi = fn_lie(bq{7},bdq{3});
fn_plot(tpsi,9);
tpsi = fn_lie(bdq{4},tpsi);
fn_plot(tpsi,10);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function fn_plot(psi,nfig)

close(figure(nfig));
figure(nfig);
hold on;

xs = 0.3;

amp = psi.amp;
base = psi.base;

left = 0;
right = max(max(base)) + 1;

vx = 1:(right-1);
vy = zeros(1,length(vx));

plot(vx, vy, 'o','markersize',6);


height = 1;

L = length(amp);

for i = 1:L
    
   a = amp(i);
   v = base(i,:);
   m = v(1);
   n = v(2);
   
   if m > n
       
       vx = linspace(m,n,100);
       vy = sqrt(((m-n)/2)^2 - (vx - (m+n)/2).^2);
       plot(vx,vy,'b-');
       text((m+n)/2 - xs,(m-n)/2,num2str(a));
       
       if abs((m-n)/2)>height
           height = abs((m-n)/2);
       end
       
   elseif m < n
       
       vx = linspace(m,n,100);
       vy = sqrt(((m-n)/2)^2 - (vx - (m+n)/2).^2);
       plot(vx,-vy,'r-');
       text((m+n)/2 - xs,(m-n)/2,num2str(a));
       
       if abs((m-n)/2)>height
           height = abs((m-n)/2);
       end
       
   elseif m == n
       
       vx = linspace(m-1/3+0.0001,m+1/3- 0.0001,100);
       vy = sqrt((1/3)^2 - (vx - m).^2);
       plot(vx, vy, 'g-');
       plot(vx, -vy, 'g-');
       text(m - xs,1/3,num2str(a));
       
   end
    
end

axis([left, right, -2*height, 2*height]);

end


function fn_print(psi)

amp = psi.amp;
base = psi.base;

L = length(amp);

if L == 0
    fprintf('0 \n');
    return;
end

if L == 1
    if amp > 0
        str = compose("%4.2f",amp);
        str = str+"("+int2str(base(1))+","+int2str(base(2))+")";
    end
else
    total_length = 0;
    str_add = compose("%4.2f", amp(1));
    str_add = str_add +"("+int2str(base(1,1))+","+int2str(base(1,2))+")";
    total_length = total_length + strlength(str_add);
    
    if total_length>50
        str_add = str_add + "\n";
        total_length = 0;
    end
    
    str = str_add;
    
    for n = 2:L
        a = amp(n);
        v = base(n,:);
        if a > 0
            str = str + "+";
        end
        
        str_add = compose("%4.2f", a);
        str_add = str_add + "("+int2str(v(1))+","+int2str(v(2))+")";
        total_length = total_length + strlength(str_add);
        if total_length>50
            str_add = str_add + "\n";
            total_length = 0;
        end
        
        str = str + str_add;
        
    end
    
end

fprintf(str);

fprintf('\n');

end


function psi = fn_bq(q,Num)

amp = zeros(Num - q,1);
base = zeros(Num - q, 2);

for i = 1:(Num - q)
    
   amp(i) = 1;
   base(i,:) = [i, i+q];
   
end

psi.amp = amp;
psi.base = base;
    
end


function psi = fn_bdq(q,Num)

amp = zeros(Num - q,1);
base = zeros(Num - q, 2);

for i = 1:(Num - q)
    
   amp(i) = 1;
   base(i,:) = [i+q, i];
   
end

psi.amp = amp;
psi.base = base;
    
end




function psi = fn_lie(psi_1, psi_2)

L1 = length(psi_1.amp);
L2 = length(psi_2.amp);

psi.amp = [];
psi.base = [];

amp1 = psi_1.amp;
amp2 = psi_2.amp;
base1 = psi_1.base;
base2 = psi_2.base;

for i1 = 1:L1
    for i2 = 1:L2
        
        vec1 = base1(i1,:);
        vec2 = base2(i2,:);
        t1 = amp1(i1);
        t2 = amp2(i2);
        tmp = fn_lie_oneterm(vec1, vec2);
        tmp.amp = t1*t2*tmp.amp;
        psi = fn_add(psi,tmp);
    end
end

end


function psi = fn_lie_oneterm(a,b)

m = a(1);
n = a(2);
p = b(1);
q = b(2);

psi0.amp = [];
psi0.base = [];

if m == n && p == q
    psi = psi0;
    return;
end

if m == n && p ~= q
    if p ~=m && q ~=m
        psi = psi0;
    else
        if p == m
            psi.amp = 1;
            psi.base = [m,q];
        else
            psi.amp = -1;
            psi.base = [p,m];
        end
    end
end

if p == q && m ~= n
    if m ~= p && m ~=q
        psi = psi0;
    else
        if m == p
            psi.amp = -1;
            psi.base = [p,n];
        else
            psi.amp = 1;
            psi.base = [m,p];
        end
    end
end

if p~=q && m~=n
    if m == q && p ~= n
        psi.amp = -1;
        psi.base = [p,n];
    elseif n == p && m ~= q
        psi.amp = 1;
        psi.base = [m,q];
    elseif m == q && n == p
        psi.amp = [1,-1];
        psi.base = [m, m; n, n];
    else
        psi = psi0;
    end
end
    
end





function psi = fn_add(psi1, psi2)

amp = [psi1.amp, psi2.amp];
base = [psi1.base; psi2.base];

psi.amp = amp;
psi.base = base;
psi = fn_combine(psi);

end


function psi = fn_combine(psi1)

amp = psi1.amp;
base = psi1.base;

Num = length(amp);

for i = 1:(Num - 1)  
    if amp(i)~=0
        vec0 = base(i,:);
        for i2 = (i + 1):Num
            vec1 = base(i2,:);
            if max(abs(vec1 - vec0)) == 0
                amp(i) = amp(i) + amp(i2);
                amp(i2) = 0;
            end
        end
    end
    
end

ind = find(abs(amp) < 1e-12);

amp(ind) = [];
base(ind,:) = [];

psi.amp = amp;
psi.base = base;

end